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ABSTRACT 

The three-dimensional (3D) coordinates of stereoscopically triangulated loops provide strong 
constraints for magnetic field models of active regions in the solar corona. Here we use 
STEREO/ A and B data from some 500 stereoscopically triangulated loops observed in four 
active regions (2007 Apr 30, May 9, May 19, Dec 11), together with SOHO/MDI line-of-sight 
magnetograms. We measure the average misalignment angle between the stereoscopic loops and 
theoretical magnetic field models, finding a mismatch of fi = 19° — 46° for a potential field 
model, which is reduced to [i = 14° — 19° for a non-potential field model parameterized by 
twist parameters. The residual error is commensurable with stereoscopic measurement errors 
i^SE ~ 8° — 12°). We developed a potential field code that deconvolves a line-of-sight mag- 
netogram into three magnetic field components (B x , B yi B z ), as well as a non-potential field 
forward-fitting code that determines the full length of twisted loops (L ss 50 — 300 Mm), the 
number of twist turns (median N twist — 0.06), the nonlinear force-free a-parameter (median 
a«4x 10~ n cm^ 1 ), and the current density (median j z w 1500 Mx cm~ 2 s _1 ). All twisted 
loops are found to be far below the critical value for kink instability, and Joule dissipation of 
their currents is found be be far below the coronal heating requirement. The algorithm devel- 
oped here, based on an analytical solution of nonlinear force-free fields that is accurate to second 
order (in the force- free parameter a), represents the first code that enables fast forward-fitting 
to photospheric magnetograms and stereoscopically triangulated loops in the solar corona. 

Subject headings: Sun: Corona — Sun: Magnetic field — Sun: EUV radiation 



1. INTRODUCTION 

This is paper IV of a series that explores the three-dimensional (3D) reconstruction of coronal loops with 
data from the two STEREO A and B spacecraft. The previous three studies focused on the 3D geometry of 
coronal loops (Aschwanden ct al. 2008b), on electron density and temperature measurements (Aschwanden et 
al. 2008c), and on active region modeling with a stereoscopic tomography method (Aschwanden et al. 2009). 
In this fourth study we concentrate on magnetic modeling of four stereoscopically observed active regions, 
which should reveal whether potential or non-potential magnetic field models fit the observed data better, 
and should provide more accurate information on magnetic fields and electric currents that are directly 
inferred from observables, and this way allow us to test theoretical magnetic field models of the solar corona. 



1 Manuscript version, 2012-Jul-4 



-2- 



The crucial benefit of this study is the inference of the true 3D magnetic field geometry, which presently can 
only be accurately measured by stereoscopic triangulation with the dual spacecraft STEREO/ A and B, and 
thus it provides a most crucial input to test theoretical magnetic field models. 

The Solar TErrestrial RElations Observatory (STEREO) was launched on October 26, 2006, consisting 
of two identical spacecraft A(head) and B(ehind) that orbit from the Earth in opposite directions around 
the Sun, with a separation angle that increases by s=s 45° per year (Kaiser 2008). Some tutorials or reviews 
on magnetic modeling using STEREO data can be found in Inhester (2006), Wiegelmann ct al. (2009), and 
Aschwanden (2011). Early methods of stereoscopic 3D magnetic field modeling have been already proposed 
in the pre-STEREO era, including loop twist measurements by fitting the projected shape of curved helical 
field lines (Portier-Fozzani et al. 2001), by fitting a linear force- free model to coronal loops geometrically 
reconstructed with solar-rotation stereoscopy (Wiegelmann and Neukirch 2002; Feng et al. 2007a), by tomo- 
graphic reconstruction constrained by the magnetohydrostatic equations (Wiegelmann and Inhester 2003; 
Ruan et al. 2008), by comparison of theoretical magnetic field models with spectro-polarimetric loop de- 
tections (Wiegelmann et al. 2005), or by using theoretical magnetic field models to resolve the stereoscopic 
correspondence ambiguity of triangulated loops (Wiegelmann and Inhester 2006; Conlon and Gallagher 2010). 
Once the STEREO mission was launched, the first stereoscopic triangulations of loops were performed from 
data of April 2007 onward, with a spacecraft separation angle of a sep > 5° (Feng et al. 2007b; Aschwanden 
et al. 2008b). First fits of linear force- free field extrapolations (using SOHO/MDI magnctograms) to stereo- 
scopically triangulated loops were found to be offset (Feng ct al. 2007b; Inhester et al. 2008). A benchmark 
test of one potential field code and 11 nonlinear force- free field codes modeling a partial active region (NOAA 
10953) observed with Hinode/SOT and STEREO/EUVI on 2007 Apr 30 revealed mean 3D misalignment 
angles of ii — 24° for the potential field code, and a range of /j, = 24° — 44° for the non-potential field 
(NLFFF) codes, which was attributed to the non-force-freeness in chromospheric heights and uncertainties 
in the boundary data (DeRosa et al. 2009). Misalignment measurements between stereoscopically triangu- 
lated loops and potential field source surface (PFSS) magnetic field extrapolations were then extended to 
four active regions with a similar finding of fi = 19° — 36° (Sandman et al. 2009). The misalignment could 
be reduced by about a factor of two for potential field models that were not constrained by observed magne- 
tograms, where the magnetic field is parameterized by buried unipolar magnetic charges (Aschwanden and 
Sandman 2010) or by a small number of 3-5 submerged dipoles (Sandman and Aschwanden 2011). These 
exercises demonstrated that potential field solutions exist that match the 3D geometry of coronal loops 
better than standard extrapolations from photospheric magnetograms. 

Modeling of non-potential magnetic fields, constrained by stereoscopically triangulated loops, however, 
is still unexplored. A first step in this direction was pioneered by Malanushenko et al. (2009), taking a 
particular model of a nonlinear force-free field with a known analytical solution (Low and Lou 1990) and 
fitting a linear force-free field individually to each field line, which yields the twist and force-free a-parameter 
for each loop individually, but not in form of a self-consistent nonlinear force-free field. In a next step, a 
self-consistent nonlinear force-free field was fitted to the en's or each loop (Malanushenko et al. 2012). In this 
study wc pursue for the first time non-potential magnetic field modeling applied to an observed set of coronal 
loops, i.e., to some 500 loops that have been stereoscopically triangulated from four different active regions. 
Since currently available NLFFF codes are very computing-intensive, which precludes an iterative fitting to 
individual loops, we developed a magnetic field parameterization that is suitable for fast forward-fitting. In 
this parameterization, the magnetic field is characterized by a superposition of point charges with vertically 
twisted fields, which is approximately force-free (to second order in a), as derived analytically (Aschwanden 
2012) and tested numerically (Aschwanden and Malanushenko 2012). Since forward-fitting requires some 
form of model parameterization, tests with physics-based models (such as twisted field components used 
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here) will be useful to explore how closely the observed loop geometry can be fitted at all, although our 
choice of parameterization represents only a subset of all possible nonlinear force-free solutions. We will test 
the validity of the divergence-freeness and force-freeness numerically and quantify it by common figures of 
merit, which can be compared with the performance of other full-fledged NLFFF codes (which, however, are 
not capable of fast forward- fitting to observations) . 

The outline of this paper is as follows: Section 2 contains the theoretical outline of a force-free field 
approximation, Section 3 presents the results of forward-fitting to some 500 stereoscopically triangulated 
loops, and Section 4 concludes with a discussion of the results. Technical details of the forward-fitting 
code are described in Appendix A and in two related documentations (Aschwanden 2012; Aschwanden and 
Malanushenko 2012). 

2. THEORETICAL MODEL AND DATA ANALYSIS METHOD 

Here we develop a novel method to determine an approximate 3D magnetic field model of a solar active 
region, using a magnetogram that contains the line-of-sight magnetic field component and a set of stereo- 
scopically triangulated loops, observed with STEREO/ A and B in extreme ultraviolet (EUV) wavelengths. 
A flow chart of this magnetic field modeling algorithm is given in Fig. 1. 

2.1. 3D Potential Field with Unipolar Magnetic Charges 

Magnetograms that contain the line-of-sight magnetic field component B z (x,y) are readily available, 
while vector magnetograph data that provide the 3D magnetic field components [B x (x, y), B y (x, y), B z (x, y)\ 
are only rarely available, are difficult to calibrate, resolving the 180° ambiguity is challenging (even with the 
recent HMI/SDO instrument), and they generally contain substantially larger data noise in the B-components 
than magnetograms. It is therefore desirable to develop a method that infers 3D magnetic field components 
from magnetogram images. Such a method was recently developed with the concept of parameterizing the 
3D magnetic field with subphotospheric unipolar magnetic charges that can be determined from an observed 
line-of-sight magnetogram, as described in Aschwanden and Sandman (2010). In the first version we neglected 
the curvature of the solar surface, because the analyzed active regions have generally a substantially smaller 
size than the solar radius. For higher accuracy of the magnetic field model, however, we generalize the 
method by including here the full 3D geometry of the curved solar surface. 

The simplest representation of a magnetic potential field that fulfills Maxwell's divergence-free condition 
(V • B = 0) is a magnetic charge that is buried below the solar surface (to avoid magnetic monopoles in the 
corona), which predicts a magnetic field B(x) that points away from the buried unipolar charge and whose 
field strength falls off with the square of the distance r, 

B(x) = 5 (^, (1) 

where B is the magnetic field strength at the solar surface directly above the buried magnetic charge, 
i"o = (#o>2/o> z o) is the subphotospheric position of the buried charge, do = yl — x\ — y^ — Zq is the depth 
of the magnetic charge, and r = [(x — xq), (y — yo), (z — zq)] is the distance vector of an arbitrary location 
x = (a;, y, z) in the solar corona (where we desire to calculate the magnetic field) from the location ro of the 
buried charge. We choose a cartesian coordinate system (x, y, z) with the origin in the Sun center and use 
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units of solar radii, with the direction of z chosen along the line-of-sight from Earth to Sun center. For a 
location near disk center (x <C 1, y <C 1), the magnetic charge depth is d ~ (1 — z ), an approximation that 
was used earlier (Aschwanden and Sandman 2010), while we use here the exact 3D distances to account for 
the curvature of the solar surface and for off-center positions of active region loops. 

Our strategy is to represent an arbitrary line-of-sight magnetogram with a superposition of N m magnetic 
charges, so that the potential field can be represented by the superposition of N m fields Bj from each magnetic 
charge j = 1, ...,N m , 

B(x) = ^B j (x)=^S j (^) 1, (2) 

3=1 3=1 

with r = [(x — Xj), (y — yj),(z — Zj)] the distance from the magnetic charge j. Since the divergence operator 
is linear, the superposition of a number of potential fields is divergence-free also, 

V-B = V-(^B,) = ^(V-B j ) = 0. (3) 

j 3 

This way we can parameterize a 3D magnetic field B (x) (Eq. 2) with 4 x N m parameters, i.e. , (xj , j/j , Zj ,Bj),j — 
l,...,N m . 

A simplified algorithm to derive the magnetic field parameters (xj,yj,Zj,Bj),j = l,...,N m from a 
line-of-sight magnetogram B z (x,y) is described in the previous study (Aschwanden and Sandman 2010). 
Essentially we start at the position (xi,yi) of the peak magnetic field B z (xi,yi) in the magnetogram image, 
fit the local magnetic field profile B x (x) (Eq. 1) to obtain the parameters of the first magnetic field component 
(xi,yi,Zi,B\), subtract the first component from the image B z (x,y), and then iterate the same procedure 
at the second-highest peak to obtain the second component (#2 , yi , £2 , Bi ) , and so forth, until we stop at 
the last component (xN m , yN m , z N m , Bw m ) above some noise threshold level. We alternate the magnetic 
polarity in the sequence of subtracted magnetic charge components, in order to minimize the absolute value 
of the residual maps. The number of necessary components typically amounts to N m s=y 100 for an active 
region, depending on the complexity and size of the active region. The geometrical details of the inversion 
of the parameters (xj,yj,Zj,Bj) of a magnetic field component from the observables (B z ,x p7 y P7 w) of an 
observed peak in a line-of-sight magnetogram at position (x p , y p ) with peak value B z and width w is derived 
in Appendix A. 



2.2. Force-free Magnetic Field Model 

With the 3D parameterization of the magnetic field described in the foregoing section we can compute 
magnetic potential fields that fulfill Maxwell's divergence-free condition and are "current-free" , 

V ' B = ° (4) 
j/c =I(VxB)=0. U 

However, magnetic potential field models do not fit observed 3D loop geometries sufficiently well (DeRosa et 
al. 2009; Sandman et al. 2009; Aschwanden and Sandman 2010; Sandman and Aschwanden 2011). Therefore, 
we turn now to non-potential magnetic field models of the type of nonlinear force-free fields (NLFFF), 

V ' B = ° (5) 
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where a(x) is a scalar function that varies in space, but is constant along a given field line, and the current 
j is co-aligned and proportional to the magnetic field B. General NLFFF solutions of Eq. (5), however, arc 
very computing-intensive and are subject to substantial uncertainties due to insufficient boundary constraints 
and non-force-free conditions in the lower chromosphere (Schrijver et al. 2006; DeRosa et al. 2009). In our 
approach to model NLFFF solutions that fit stereoscopically triangulated loops we thus choose approximate 
analytical NLFFF solutions that can be forward-fitted to the observed 3D loop geometries much faster. 
The theoretical formulation of the force-free magnetic field approximation suitable for fast forward-fitting is 
described in Aschwanden (2012) and numerical tests are given in Aschwanden and Malanushenko (2012). 

The essential idea is that we generalize the concept of buried point charges we used for the representation 
of a potential field, by adding an azimuthal magnetic field component B v (p) that describes a twist around 
a vertical axis. The concept is visualized in Fig. 2. An untwisted flux tube can be represented by parallel 
field lines B s (p) along an axis s with variable distances p from a given axis s (Fig. 2, left). A uniformly 
twisted flux tube with some azimuthal component B ip (p) around the twist axis s (Fig. 2, middle) has the 
following analytical solution for a force-free field (Gold and Hoyle 1958; Priest 1982, p. 125; Sturrock 1994, 
p.216; Boyd and Sanderson 2003, p.102; Aschwanden 2004, p.216), 

BM = • (•) 

where 6 is a parameter that quantifies the number N tW i S t of full twist turns over a (loop) length L, 

2irN twist 

b =^—> (8) 
and is related to the force-free parameter a in the force-freeness condition (V x B) = a(x)B by, 

a{p) = (TT&V) ■ (9) 

Thus, we see that a has a finite value a(p = 0) = 26 at the twist axis and drops monotonically with increasing 
distance p from the twist axis. 

In a next step we generalize the concept of a cylindrical twisted flux tube (Fig. 2, middle) to a point 
charge with an associated vertically twisted field (Fig. 2, right). The main difference to the twisted flux 
tube, which has a constant cross-section and thus a constant magnetic flux along the cylinder axis, is the 
quadratically decreasing field strength of the longitudinal field component B s (s) along the twist axis with 
distance s, to conserve the magnetic flux. It can be shown analytically, that this topology has the following 
approximate force-free and divergence-free solution (satisfying Eqs. 4 and 5), if we neglect second-order and 
higher-order terms of [6rsin(6*)] (Aschwanden 2012), 

B r (r, 6) = B (*) 1 , , (10) 

\r 2 J (1 + 6 2 r 2 sin 2 6) 

, f\ \ I d 2 \ br sin 9 , . 

s e M)~o, (12) 
. 26 cos ,„„. 

(1 + b 2 r 2 sin 0) 
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Since b oc a for 9 <C 1 (Eq. 13), we refer to this approximation as second-order accuracy in a for short. 
Thus we can characterize a point charge with a vertically twisted field with 5 parameters: {Bj, Xj, yj, Zj, 
where the force-free parameter ctj = a(p = 0) at the twist axis is related to bj by a., = 2b j. 

In analogy to the superposition principle we used to construct an arbitrary potential field based on a 
number of N m magnetic charges, we apply the same superposition rule to a number of N m point charges 
with vertical twist, 

N m 

B(x) = ^B j (x), (14) 
i=i 

and we can parameterize an arbitrarily complex non-potential field with 5 x N m free parameters. The 
superimposed field B, which consists of a linear combination of quasi- force- free field components (to second 
order in a) is not exactly force-free, but it can be shown that the NLFFF approximation expressed in 
Eq. (10-13) is divergence-free to second-order accuracy in the term [br sm(9)], and force-free to third-order 
accuracy (Aschwanden 2012), which we call "quasi-force-free" here. The superposition of multiple twisted 
sources, which have each a twist axis with a different location and orientation (due to the curvature of the 
solar surface) , needs to be transformed into a common cartesian coordinate system, which is derived in detail 
in Section 2.4 in Aschwanden (2012). 

We show a few examples of this non-potential magnetic field model in Fig. 3, for a unipolar, a dipolar, 
and a quadrupolar configuration. More cases are simulated in Aschwanden (2012). Three cases are shown 
for the non-potential case with twist (Fig. 3, right panels), which degenerate to the potential field case if the 
twist parameter is set to zero (otj = 0) (Fig. 3, left panels). 

In our data analysis we perform a global fit of our non-potential field model to a set of stereoscopically 
triangulated loops, Typically we use a magnetic field model with N m < 200 magnetic field components 
(which has < 800 potential field parameters that can directly be inverted from the line-of-sight magnetogram 
B z (x,y 7 z p h), and < 200 non-potential field parameters fx, , which need to be forward-fitted to the data, 
consisting of the stereoscopic 3D coordinates of some Nl ~ 70 — 200 loops per active region. Each forward- 
fit minimizes the average misalignment angle (along and among the loops) between the theoretical model (i.e., 
twisted fields) and the observed loop (with stereoscopically triangulated 3D coordinates). The misalignment 
angle at a position x is defined as (Sandman et al. 2009; Aschwanden and Sandman 2010; Sandman and 
Aschwanden 2011): 

^-■""■(I^m)' (15 » 

We derive a characteristic misalignment angle from measurements at 10 equi-spaced postions x along each 
observed loop segment, where a unique theoretical field line B t /j eo (x) is calculated at each loop position. 
Generally, this yields 10 theoretical field lines that intersect a single loop at the 10 chosen positions, which 
degenerate to one single theoretical field line in the case of a perfectly matching model. The range of 
misalignment angles ignores the 180°-ambiguity and is defined in the range of \i = [0°, 90°]. Technical details 
of the numerical algorithm of our forward-fitting code to the stereoscopic loop coordinates are provided in 
Aschwanden and Malanushcnko (2012). 

2.3. Stereoscopical Triangulation of Loops 

In order to obtain the 3D coordinates of coronal loops we require one (or multiple wavelength filter) image 
pairs recorded with the two spacecraft STEREO/ A(head) and B(ehind), at some spacedraft separation angle 
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that is suitable for stereoscopy. Here we will use STEREO data from the first year (2007) of the mission, 
when the spacecraft separation was a sep < 45°. An example of such a stereoscopic pair of EUVI images 
recorded at a wavelength of 171 A is shown in Fig. 1 (top right). Loop recognition requires a segmentation 
algorithm that traces curvi-linear features. There are visual/manual methods as well as automated methods, 
which have been compared in a benchmark study (Aschwanden et al. 2008a). Although automated detection 
algorithms have improved to the level or approaching visual perception for single images (Aschwanden 2010), 
an automated algorithm for dual detection of corresponding loop patterns in stereoscopic image pairs has 
not yet been developed yet, so that we have to rely on manual/visual detection. The 3D coordinates of 
such visually detected loops have been documented in previous publications (Aschwanden 2008b, c, 2009; 
Sandman et al. 2009; Aschwanden and Sandman 2010; Sandman and Aschwanden 2011). 

The technique of loop detection falls into the category of image segmentation algorithms. The easiest 
way to segment coronal loops in an EUV image is to apply first a highpass filter (Fig. 1), which enhances 
fine structures, in particular coronal loops with narrow widths. Stereoscopic triangulation can then be ac- 
complished by (1) co- aligning a stereoscopic image pair into epipolar coordinates (where the stereoscopic 
parallax effect has the same direction in both images), (2) tracing a loop visually in a STEREO/ A image, 
(3) identifying the corresponding counterpart in the STEREO/B image using the expected projected alti- 
tude range; (4) stereoscopic triangulation of the loop coordinates (sl[s],|/£[s],^l[s]) (along the loop length 
coordinate s), and (5) coordinate transformations into an Earth-Sun heliographic coordinate system. Here 
we have to transform STEREO/ A coordinates into the SoHO/MDI coordinate system, in order to enable 
magnetic modeling with MDI magnetograms. The height range of stereoscopically triangulated loops does 
generally not exceed 0.15 solar radii, due to the decrease of dynamic range in flux for altitudes in excess 
of one hydrostatic scale height. Here we used the loop tracings and stereoscopic measurements from three 
different EUVI channels combined (171, 195, 284 A). The method of stereoscopic triangulation is described 
in detail in previous papers (Aschwanden 2008b, c, 2009). 

What is also important to know when we attempt magnetic modeling of stereoscopically triangulated 
coronal loops is the related uncertainty. The parallax effect occurs mostly in east- west direction, which implies 
that loop segments with North-South direction can be most accurately triangulated, while the accuracy 
decreases progressively for segments with longer components in east-west direction (Aschwanden 2008b). 
Another source of error is mis-identification or confusion of dual loop counterparts in stereoscopic image 
pairs, an uncertainty that can only be quantified in an empirical manner. Such empirical estimates of 
stereoscopic errors were evaluated by measuring deviations from parallelity of loop pairs triangulated in 
near proximity, which yielded an estimated uncertainty in the absolute direction of use — 8° — 12° for four 
analyzed active regions (Aschwanden and Sandman 2010). So, even when we achieve a theoretical magnetic 
field model that perfectly matches the observed loops, we expect a residual misalignment error of /i « 10° 
due to the uncertainty of stereoscopic triangulation. 

2.4. 3D Magnetic Field Representation 

The input from the line-of-sight magnetogram provides the longitudinal component B z (x, y, z p h) at the 
photospheric level r = 1, which corresponds to the z-coordinate z p h = \J (\ — x 1 — y 2 ). The decomposition 
of the magnetogram data into unipolar magnetic charges (Section 2.1 and Appendix A) yields the full 3D 
vector field, B p (x) = [B x (x, y, z p h) 7 B y (x, y, z p h), B z (x, y 7 z p h)} on the solar surface at the photospheric level 
z^, which is a potential field by definition, since the derived field consists of a superposition of potential 
fields from single magnetic charges (Eq. 2). 
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The nonlinear force-free magnetic field model (Section 2.2), which we fit to the data by optimizing the 
misalignment angles with our chosen parameterization, defines a non-potential field B(x) = [B x (x, y, z), 
B y (x, y, z), B z (x, y, z)], which can be extrapolated at any point in the corona by stepping along a field line 
[x(s), y(s), z(s)] in small incremental steps As, where the cosines of the cartesian magnetic field components 
(B x , B y ,B z ) define also the direction of the field line, so that it can be iteratively calculated with 

x(s + As) = x(s) + As[B x (s)/B(s)]p 

y(s + As) =y(s) + As[B y (s)/B(s)]p . (16) 
z(s + As) = z(s) + As[B z (s)/B(s)]p 

where p — ±1 represents the sign or polarization of the magnetic charge. 

The force-free parameter a(x,y, z) of our forward-fitted field B(x,y,z) at any given point of space in 
the computational box (xi, yj, Zk) can be numerically computed for each of the three vector components of 
(V x B), 

1 (V x B) x 1 (8B Z 8B y \ 
«-( x ) = — 5 = TIB" I — ~ ) > ( 17 ) 

(18) 



1 (V x B) x 
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a ^ = -,-^B— = 4^1 {-Ox- -di) ' (19) 

using a second-order scheme for the spatial derivatives, i.e., dB x /dy — (-Bi+ij.fc — Bi-ij t k)/2(yi+i — J/i-i)- 
In principle, the three values a Xl a y , a z should be identical, but the numerical accuracy using a second- 
order differentiation scheme is most handicapped for those loop segments with the smallest values of the 
i3-component (appearing in the denominator), for instance in the a z component oc (l/B z ) near the loop 
tops (where B z ss 0). It is therefore most advantageous to use all three parameters a x , a y , and a z in a 
weighted mean, 

a x w x + a y w y + a z w z 
a = — , (20) 

w x + w y + w z 

but weight them by the magnitude of the (squared) magnetic field strength in each component, 

w x = Bl , w v = B 2 y , w z = B 2 Z , (21) 

so that those segments have no weight where the _B-component approaches zero. This numerical method 
was found to render the a-parameters most accurately (Aschwanden and Malanushenko 2012), while the 
analytical approximation (Eq. 13) breaks down near loop tops (where B z w 0). 

The computation of the current densities j(x, y, z) follows then directly from Eq. (5), 

j(x,y,z)/c = a(x,y,z)B(x,y,z) . (22) 

For instance, a photospheric (vertical) current map, i.e., j z (x,y, z p h), as shown in in Fig. 1 (bottom panel), 
can be calculated with Eq. (22) for the photospheric level at z p h — — x 1 — y 2 ). 



3. OBSERVATIONS AND RESULTS 

We present now the results from four active regions, which all have been modeled with different magnetic 
field models before also (DeRosa et al. 2009; Sandman et al. 2009; Aschwanden and Sandman 2010; Sandman 
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and Aschwandcn 2011) and can now be compared with our modeling results. A summary of the four active 
regions, including the observing time, heliographic position, STEREO spacecraft separation, number of 
stereoscopically triangulated loops, maximum magnetic field strength, and magnetic flux is given in Table 
1, which is reproduced from Aschwanden and Sandman (2010). 



3.1. 3D Potential Magnetic Field 



The first step of our analysis is the decomposition of the observed SOHO/MT)! magnctograms of four 
active regions into N m magnetic charges (Eq. 2), which yields the 4 x N m parameters [Bj, Xj, yj, Zj], j = 
1, N m for the 3D parameterization B(x) of the magnetic field. In the remainder of the paper we will label 
the four active regions observed at the given times as follows (see also Table 1): 



(A) for active region 10953, 2007-Apr-30, 23:00 UT, 

(B) for active region 10955, 2007-May-9, 20:30 UT, 

(C) for active region 10953, 2007-May-19, 12:40 UT, 

(D) for active region 10978, 2007-Dec-ll, 16:30 UT. 



The pixel size of MDI magnetograms is 2" and the MDI magnctograms have been decomposed into Gaussian- 
like peaks (each one corresponding to a buried unipolar magnetic charge) . The chosen field-of- views (FOV) 
of the four active regions are (in units of solar radii): 



FOV{A) = [ Xl = 
FOV(B) = [ Xl = 
FOV(C) = [xi = 
FOV(D) = [ Xl = 



0.65, x 2 = +0.05,yi 
0.55, x 2 = -0.25, yi 
0.22, x 2 = +0.18, yi 
0.28, x 2 = +0.12, yi 



-0.50,y2 = +0.20], 
-0.25,y2 = +0.05], 
-0.15,2/2 = +0.25], 
-0.35,2/2 = +0.05]. 



So, all four active regions are located near disk center [x = 0, y = 0], but some extend out to 0.65 solar radii 
(FOV A). We show the observed magnetograms B z (x,y) for the 4 active regions in Fig. 4 (left column), 
the corresponding model maps B z (x,y) built from N m — 200 magnetic charge components (Fig. 4, middle 
column), and the difference between the observed and model maps (Fig. 4, right column), all on the same grey 
scale, so that the fidelity of the model can be judged. The residual fields of the decomposed magnetograms 
after subtraction of N m — 200 components are 



B z /B max = -0.0003+0.0097 for (A), 
B z /B max = -0.0011 + 0.0074 for (B), 
B z /B max = -0.0004+0.0164 for (C), 
B z /B max = -0.0008 + 0.0189 for (D), 



(with the maximum field strengths B max listed in Table 1). Therefore, our decomposition algorithm (Ap- 
pendix A) represents the observed magnetic field down to a level of w 1% — 2% residuals. Note that this 
method yields the 3D components of the magnetic field, [B x ,B y , B z ], with the accuracy specified for the case 
of a potential field, while the accuracy for a non-potential field components (i.e., horizontal components B y 
and B z ) probably does not exceed a factor of two of this accuracy (i.e., w 2% — 4%), which is still better 
than the accuracy of currently available vector magnetograph data (in the order of w 10% for strong fields 
and worse for weak fields; Marc DeRosa, private communication). 
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3.2. Forward-Fitting of 3D Force-Free Field 

In Fig. 5 we show the comparison of Nl = 200 stereoscopically triangulated coronal loops for active 
region A (2007 Apr 30), obtained from STEREO spacecraft EUVI/A and EUVI/B data, with the potential- 
field model composed from N m — 100 buried magnetic charges (Eq. 2). For each modeled magnetic field 
line we have chosen the midpoint of the observed STEREO loop segment as the starting point, from which 
we extrapolate the potential field in both directions for segments (red curves) that are equally long as the 
observed loop segments (blue curves). A histogram of the median misalignment angles is shown in Fig. 5 
(bottom), measured at 10 positions along each loop segment, which has an average (A) of = 39.9°, a 
median (M) of fi = 29.7°, or a Gaussian peak (P) fit with a centroid and width of \i = 29.1 ± 21.6°. So, 
there is broad distribution of misalignment angles in the range of fi rts 10° — 60°, which clearly indicates that 
a potential field model is not a very good fit for these observables. 

In Figs. 6-9 we present the main results of our force-free field forward-fitting code applied to the stereo- 
scopically triangulated coronal loops data, for all four active regions A, B, C, and D. The 3D coordinates of 
the observed stereoscopic loops (Figs. 6-9; blue curves) and forward-fitted theoretical field lines (Figs. 6-9: 
red curves) are shown in three projections, in the x — y plane (in direction of the line-of-sight), as well as in 
the orthogonal x — z and z — y planes. The theoretical field lines are extrapolated to equally long segments as 
the observed EUV loops. The field lines have been computed with a step size of As = 0.002 solar radii (i.e., 
w 1500 km), which is identical to the spatial resolution of the magnetogram from SOHO/MDI (2.0" w 1500 
km), and commensurable with the spatial resolution of STEREO/EUVI (i.e., 2 pixels with a size of 1.6", 
which is w 2000 km). 

Fig. 5 and Fig. 6 show a direct comparison of modeling potential and non-potential fields for active 
region A, which results into a broad Gaussian distribution of fj, = 29.1° ± 21.6° for the potential field, and 
to a much narrower distribution of fj, = 16.7° ± 6.7° for the non-potential field model. A online-movie 
that visualizes the forward-fitting process from the initial guess of a potential field model to the best-fit 
non-potential model is also included in the online electronic supplementary material of this paper. 

3.3. Misalignment Statistics 

A quantitative measure of the agreement between a theoretical magnetic field model and the 3D coor- 
dinates of observed coronal loops is the statistics of misalignment angles, which is shown for all four active 
regions in the bottom panels of Figs. 6-9. We find the following Gaussian distributions of misalignment 
angles for each of the four active regions" 

H = 16.7° ± 6.7° for A (2007 April 30; Fig. 6), 
/i = 16.0° ± 6.0° for B (2007 May 9; Fig. 7), 
H = 19.2° ± 6.7° for C (2007 May 19; Fig. 8), 
fj, = 14.3° ± 5.6° for D (2007 Dec 11; Fig. 9), 

so they have a most frequent value of /j, w 15° — 19°, which matches very closely the estimated uncertainty of 
stereoscopic errors (based on the parallelity of loops in near proximity), which was found in a similar range 
of Use ~ 7° — 10° (Aschwanden and Sandman 2010). Online-movies that visualizes the forward-fitting 
process from the initial potential field model to the best-fit non-potential model are also included in the 
online electronic supplementary material of this paper. 
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A comparison of our results with previous magnetic modeling of the same four active regions is compiled 
in Table 2. The improved misalignment statistics obtained with our code in the range of fj, w 14° — 19° 
is about two times smaller than what was obtained with earlier NLFFF models (n[NLFFF] w 24° — 44°; 
DeRosa et al. 2009), and also about two times smaller than potential field models, i.e., with potential source 
surface models (PFSS: /i s=s 19° — 36°, Sandman et al. 2009), with potential field models with unipolar charges 
(/x w 16° — 26°, Aschwanden and Sandman 2009), also measured with our new code described in this paper 
ss 18° —42°). Thus, we find that our non-potential (force-free) field model clearly yields a better matching 
model than any previous potential or non-potential field model. 

The mean value of the misalignment angle \i in an active region depends also on the number N m 
of magnetic charges that have been used in the model, which is shown in Fig. 10 (right panels). We 
repeated the forward-fitting for N m = 1,2,4,10,20,50,100,200 magnetic charges and find that the mean 
misalignment angle fj,(N m ) generally improves (or decreases) with the number of magnetic components, up 
to N m sa 100, while for larger numbers (i.e., N m = 200) a diminuishing effect sets in (probably due to less 
efficient convergence in forward- fitting with a too large number of variables). We show the dependence of the 
mean misalignment angle fi(N m ) for both the potential field model (P) and the force-free field non-potential 
model (NP) in Fig. 10 (right panels), which reveal interesting characteristics how potential-like an active 
region is. Clearly, active region B (2007 May 9) is the most potential one, while active region C (2007 May 
19) is the most non-potential one, where a potential field model docs not fit at all, regardless how many 
magnetic components are used. A similar result was also obtained in the study of Aschwanden and Sandman 
(2010), where the potential-like active region B was associated with the lowest GOES-class level (A7), while 
the most non-potential-likc region C was exhibiting a GOES-class CO flare. 



3.4. Force-Free a and Electric Current j z Maps 

Maps of the magnetic field components B x (x,y), B y {x,y), and B z (x,y) obtained from the decomposition 
of line-of-sight magnetograms into 100 buried magnetic charges are shown for one (i.e., region A of 2007 Apr 
30) of the four analyzed active regions in Fig. 11. Note that only the line-of-sight component map B z (x,y) 
is directly observed (with SOHO/MDI), as shown in the bottom right panel in Fig. 11, while the other two 
component maps B x (x,y) and B y (x,y) are inferred from the spherical symmetry of the magnetic field of 
point charges assumed for a potential field (Eq. 1). These inferred field components define our potential field 
solution B p (x, y). 

Our non-potential field model adds additional azimuthal field components B v {x,y) that are only con- 
strained by fitting the stereoscopically triangulated loops. We show the difference between the non-potential 
and potential field components for the photospheric level in the middle column of Fig. 11. We notice that 
the difference between the non-potential and potential solution is essentially contained in the transverse field 
components B x (x, y) and B y (x, y), because the line-of-sight component B z (x, y) is an observed quantity and 
has to be matched by any model. The difference maps shown in Fig. 1 1 confirm that our model retrieves the 
non-potential field mostly from the transverse field components that arc not measured with a line-of-sight 
magnctogram. 

From our model we can also derive a photospheric force-free a(x, y) map (top right panel in Fig. 11) and 
a photospheric electric current density map j z (x,y) (middle right panel in Fig. 11), based on Eq. 22. The 
magnetograms and inferred current maps are shown for active region (A) in Fig. 11, and for the other active 
regions (B), (C), and (D) in Fig. 12. These inferred maps reveal that the locations of significant electric 
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currents are located in strong field regions, and exhibit as much small-scale fine structure as the magnetic 
field line-of-sight component B z (x,y). 

Statistics of obtained parameters of our non-potential field model is shown in form of histogrammed 
distributions in Fig. 13, lumped together from all 454 stereoscopically triangulated loops of our four analyzed 
active regions. The histograms are shown in lin-lin (Fig. 13 left) as well as in log-log representation (Fig. 13 
right) and can be characterized by powerlaw distributions. The median values are: L = 163 Mm for the 
length of the extrapolated field lines that were extrapolated through the midpoint of the observed loops, 
Ntwist — 0.06 for the number of twisted turns along the extrapolated full loop lengths; a = 4 x 1CT 11 cm -1 
for the force-free a-parameter, and \j z \ = 1500 Mx cm~ 2 s _1 for the electric current density. The statistics 
of misalignment angles is given in Fig. 10 and Table 2. 



3.5. Divergence-Freeness and Force-Freeness 

In order to test the accuracy of our analytical force-free model, which represents an approximation to a 
truly force- free field with an accuracy to second order (in a), it is useful to calculate some figures of merit. 
The divergence-freeness V • B = can be compared with the field gradient Bj Ax over a pixel length Ax, 



I f l(V-B)| 2 ir 



Similarly, the force-freeness can be quantified by the ratio of the Lorentz force, (j x B) = (V x B) x B to 
the normalization constant B 2 /Ax, 

1 f |(VxB)xBp 
Lf -Vj v \&/*x? (M) 

where B — |B|. These quantities, integrated over a computational box that covers the field-of- views of an 
active region and extends over a height ragne of Ah = 0.15 solar radii, were found to be for the four active 
regions A, B, C, and D: 



(A) L d = 

(B) L d = 

(C) L d = 

(D) L d = 



0.5 x 10 
1.1 x 10 
1.4 x 10 
1.0 x 10 



4 and Lf 
4 and Lf 
1 and Lf 
4 and Lf 



5 x 10~ 4 
7 x 10~ 4 
1.8 x lO" 1 
1.0 x 10~ 4 



In comparison, figure of L d w (8 ± 5) x 10 4 and Lf w (24 ± 23) x 10 4 were found for simulated cases 
(Aschwanden and Malanushenko 2012), and L d w (0.2 ± 8) x 10~ 4 and L f w (6 ± 2) x 10~ 4 for forward- 
fitting to the Low and Lou (1990) model, which is an analytical exact force-free solution (Aschwanden and 
Malanushenko 2012). Our values do not exceed figures of merits quoted from other NLFFF codes, i.e., 
L d w 0.07 and I/«0.1 (Schrijver et al. 2006; Table III therein). For sake of convenience we evaluated these 
figures of merit in a rectangular box tangential to the solar disk at disk center, extending over a height range 
of h = 1.00 — 1.15, which covers the strong field regions only for active regions near disk center. Of course, 
the widely used definition of divergence-freeness (Eq. 23) and force-freeness (Eq. 24) is proportional to the 
square of the normalization length scale L, i.e., L d oc L 2 and Lf oc L 2 , and thus would yield larger values 
for typical loop lengths or box sizes L >• Ax. 
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4. DISCUSSION 

4.1. The Potential Field Model with Unipolar Magnetic Charges 

As a first step we derived a potential field model of an active region. The knowledge of the potential 
field is a useful starting point for reconstructing a non-potential field model. In the asymptotic limit of small 
currents, i.e., for small values of the force- free parameter (a H 0), the force- free field solution automatically 
converges to the potential field solution. Furthermore, since our analytical force-free field approximation is 
accurate to second order (in a), the highest accuracy is warranted for a small non-potentiality, say up to 
one full turn of twist along a loop length. The fact that we measured small twists with a statistical median 
of Nt W i S t ~ 0.06 for the four analyzed active regions, justifies the neglect of second- and higher-order terms 
in our analytical approximation. 

Generally, vector magnetograph data are required to uniquely define the complete 3D magnetic field 
B(x) at the photospheric boundary. In constrast, potential-field codes can extrapolate a unique 3D field 
solution from a single B-componcnt only, i.e., from the line-of-sight component B z (x). Consequently, in 
our method of superimposed fields from buried magnetic charges, the potential field solution is unique in 
every point of space, as well as at the boundaries (neglecting exterior magnetic charges). The potential- 
field solution fixes also the line-of-sight component of the non-potential solution, B z p = B P , because this 
component is a direct observable that must be matched with every model, and thus is identical for both the 
potential and the non-potential model, while the transverse components represent the only free parameters 
for a non-potential solution, i.e., B^ p ^ B p and By P ^ B p . 

In this study we derived an algorithm that accurately calculates a potential magnetic field by properly 
including the curvature of the solar surface and for locations away from the solar disk center. We are not 
aware that such a potential-field code, defined in terms of buried unipolar magnetic charges and deconvolved 
from an observed magnetogram (Aschwanden and Sandman 2010), has been developed elsewhere, although 
similar parameterizations have been used in magnetic charge topology (MCT) models (e.g., Longcope 2005). 
It would be interesting to compare its performance with other potential field codes, such as with the Green's 
function method (Sakurai 1982), the eigenfunction expansion method (Altschulcr and Newkirk 1968), or the 
potential field source surface (PFSS) code (e.g., Luhmann et al. 1998). Our magnetic charge decomposition 
method is related to the Green's function method, but differs in the discretization of discrete magnetic 
elements of various strengths and variable location, while the Green's method uses a regular surface grid. 
Magnetic elements with sub-pixel size represent no particular problem (except that no more than one element 
can be resolved per pixel), because our method fits their radial field above the photosphere, but treats them 
as point charges below the photosphere. However, the resolution of our method is somehow limited by 
the depth of the point charges (which introduces finite numerical errors in the calcuation of V • B and 
V x B), as well as by the number of point charges (which represent a sensitivity threshold for weak magnetic 
field sources). Thus, smaller grid pixel sizes and larger number of magnetic field charges can enhance the 
accuracy of the solutions, but are more demanding regarding computation times. Similarly, the resolution 
of the Green's function method is limited by the grid size, which is either given by the measurements, or 
limited by computation ressources. 
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4.2. Non-Potential Fields with Twisted Loops 

Models of twisted flux tubes have been applied abundantly in solar physics, e.g., to braiding of coronal 
loops (Berger 1991), to prominences (Priest et al. 1989, 1996), to sigmoid-shaped filaments (Rust and Kumar 
1996; Pevtsov et al. 1997), to emerging current-carrying flux tubes (Leka et al. 1996; Longcope and Welsch 
2000), or to turbulent coronal heating (Inverarity and Priest 1995). Magnetic structures are believed to 
be twisted and current-carrying before they emerge at the solar surface (Leka et al. 1996), many active 
region loops are observed to have a visible twist (see Section 6.2.4 in Aschwanden 2004), and filaments or 
prominences become unstable and erupt once their twist exceeds a critical angle of a few full turns due 
to the kink instability or torus instability (Fan and Gibson 2003, 2004; Torok and Kliem 2003; Kliem et 
al. 2004) . The twist of a magnetic structure is therefore an important indicator for its stability or transition 
to an instability, followed by the dynamic evolution that leads to eruptive flares and coronal mass ejections 
(CME). Our method allows us to determine the number of twisting turns of most stereoscopically triangulated 
coronal loops directly, and thus provides a reliable diagnostics on its stability and the amount of electric 
current the loop carries. The statistics in Fig. 13 shows that the number of twisting turns has a range up to 
Ntwist S 0.25 (Fig. 13), which is far below the critical limit of N tW ist ~ 1-2 — 2.4 full turns required for the 
onset of the kink instability in a force- free magnetic field (Mikic et al. 1990; Torok and Kliem 2003). 

The amount of twisting clearly varies among different active regions. The improvement in the misalign- 
ment angle between a potential and a non-potential model characterizes the potentiality and free energy 
of an active region. The most potential active region is B (2007 May 9), where we achieve only a slight 
improvement of A/i = /i p — fi NP = 19° — 16° = 3° (Table 2, for the Gaussian peak of the distribution). Also 
active region D (2007-Dec-ll) is moderately potential, where we find A/x = /i P — [i NP = 21° — 15° = 6°. 
Significant non-potentiality is found for active region A (2007-Apr-30), where the improvement amounts to 
A/i = fi p - ll np = 29° - 17° = 12°. The strongest non-potentiality is found for active region C (2007-May- 
19), where we achieve an improvement of A/i = ji p — ji NP = 46° — 20° = 26°. Therefore, the misalignments 
reduce by an amount of A/i = 3°..., 26° for these four active regions. During the time of observations, active 
region C indeed featured a GOES-class CO flare, which explains its non-potentiality. The remaining amount 
of misalignment, in the order of /i = 15° — 20° is attributed partly to stereoscopic measurement errors, which 
were estimated to use ~ 8° — 12°, and partly to our particular parameterization of a force-free field model, 
which is optimally designed for vertically twisted structures. 



4.3. Force-Free a Parameter and Electric Currents 

For the force-free a parameter we find a distribution extended over the range of a < 20 x 10" 11 cm -1 
(Fig. 13), with a median of a w 4 x 10" 11 cm -1 . Using a sheared arcade model, where the force-free a 
parameter is defined as a — Ltan/i (with L the length of the loops), a range of a — (0.6 — 13) x 10" 11 
cm -1 was found for the same dataset, based on the residual misalignment of /i « 7° — 13° attributed to non- 
potentiality (Table 4 in Sandman and Aschwanden 2011). Modeling of stereoscopically triangulated loops 
with a linear force-free model yielded also a similar range of a ss (2 — 8) x 10" 11 cm -1 (Feng et al. 2007b). 
Thus, the range of our determined a parameters agrees well with these three studies. 

The current densities have been determined in a range of j ~ 10 " 2 — 10 Mx cm" 2 s , with a median 
of j s=y 1500 Mx cm" 2 s" 1 (Fig. 13), for spatial locations with B > 100 G at the photospheric level. For 
comparison, Leka et al. (1996) determined similar currents of j z « (21 — 75) Mx cm" 2 s" 1 (i.e., 7—25 
mA m" 2 , see Table 3 in Leka et al. 1996), using vector magnetograph data at the Mees Solar Observatory 
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measured in emerging bipoles in active regions. In principle, the knowledge of the magnitude of the electric 
current density j allows one to estimate the amount of Joule dissipation, 

E H = — , (25) 
a 

where a s=s 6 x 10 16 s _1 is the classical conductivity for a T k> 2 MK hot corona. Our measured maximum 
current density of \j z \ < 10 4 Mx cm -2 s _1 yields then a value of E H < 2 x 10 -9 erg cm -3 s _1 for the 
volumetric heating rate. The corresponding Poynting flux Fjj for a field line with a density scale height of 
A ss 10 10 cm (corresponding to the thermal scale height at a temperature of T = 2 MK), is 

F H = E H A(T) « li) w E H (-^) , (26) 

which yields Fh < 20 erg cm -2 s _1 , which is far below the heating requirement for active regions (Fh ~ 
10 5 — 10 7 erg cm -2 s _1 ; Withbroe and Noyes 1977). Thus, Joule dissipation is far insufficient to heat coronal 
loops in active regions, according to our current measurements. Obviously, more energetic processes, such as 
magnetic reconnection with anomalous resistivity in excess of classical conductivity is needed. Nevertheless, 
our technique of measuring current densities from the twist of observed loops may provide a useful diagnostic 
where currents are generated in active regions, e.g., near neutral lines with large gradients in the magnetic 
field, or in areas with large photospheric shear that produces magnetic stressing. 

4.4. Benefits of Non-Potential Magnetic Field Modeling 

Modeling of the magnetic field with a non-potential model matches the observed loop geometries about 
a factor of two better than potential field models, as demonstrated in this study (Table 2). While our 
attempt of forward-fitting of twisted fields to stereoscopically triangulated loops represents only a first 
step in this difficult problem of data-constrained non-potential field modeling, we expect that future non- 
potential models with similar parameterizations will be developed, which yield the 3D magnetic field with 
high accuracy, satisfying Maxwell's equations of divergence-freeness and force-freeness. There are a number 
of benefits that will result from the knowledge of a realistic 3D magnetic field, of which we just mention a 
few: 

1. The spatial distribution of coronal DC current heating can be inferred from localizations of non- 
potential field lines with a high degree of twist and current density (e.g., Bobra et al. 2008; DeRosa et 
al. 2009; Su et al. 2011). 

2. The free (non-potential) energy Ef ree = Em — Ep can be determined from the difference of the non- 
potential (Em) and potential (Ep) energy in a loop. 

3. Correlations between the free energy Ef ree and the soft X-ray brightness Isxr can quantify the vol- 
umetric heating rate Eh and pointing flux Fjj in loops and active regions, and be related to the flare 
productivity (e.g., Jing et al. 2010; Aschwanden and Sandman 2010). 

4. The helicity and its evolution can be traced in twisted coronal loops, which can be used to diagnose 
whether helicity injection from below the photosphere takes place (e.g., Malanushenko 2011b). 

5. Hydrodynamic modeling of coronal loops, which is mostly done in a one-dimensional coordinate (s) 
along the loop, requires the knowledge of the 3D geometry [x(s),y(s), z(s)} 7 which otherwise can only 



- 16 - 



be obtained from stereoscopy. The inclination of the loop plane determines the ratio of the observed 
scale height to the effective thermal scale height, and thus the inferred temperature profile T e (s). Also 
the inference of the electron density n e (s) along the loop depends on the column depth of the line- 
of-sight integration (e.g., Aschwanden et al. 1999). Most loop scaling laws depend explicitly on the 
loop length L (e.g., the Rosner- Tucker- Vaiana law, Rosner et al. 1978), which can only be reliably 
determined from magnetic models. 

6. Electron time-of-flight measurements require the trajectory of electrons that stream along field lines 
from coronal acceleration sources to the photospheric footpoints. The location of the coronal accel- 
eration can be determined from the observed energy-dependent time delays in hard X-rays and the 
knowledge of the pitch angle of the electron and the magnetic twist of the field line, which can only be 
obtained from magnetic models (e.g., Aschwanden et al. 1996). 

7. Modeling of entire active regions require models of the 3D magnetic field, the temperature T e (s), and 
density n e (s), using hydrostatic or hydrodynamic loop models, which can reveal scaling laws between 
the volumetric heating rate E H , hydrodynamic loop parameters (n e ,T e ,p, L), and magnetic parameters 
(B), (e.g., Schrijver et al. 2004; Warren & Winebarger 2006, 2007; Warren et al. 2010; Lundquist et 
al. 2008a,b). 

8. Nonlinear force- free modeling of active regions and global coronal fields can establish better lower 
boundary conditions for modeling of the heliospheric field (e.g., Petrie et al. 2011). 

5. CONCLUSIONS 

There are two classes of coronal magnetic field models, potential fields and non-potential fields, which 
both so far do not fit the observed 3D geometry of coronal loops well. Nonlinear force-free fields (NLFFF) 
are found to match the data generally better than linear force- free fields (LFFF). The 3D geometry of 
coronal loops is most reliably determined by stereoscopic triangulation, as it is now available from the twin 
STEREO/ A and B spacecraft, and we made use of a sample of some 500 loops observed in four active 
regions during the first year of the STEREO mission (with spacecraft separation angles in the range of 
&sep ~ 6° — 43°). We forward-fitted a force- free approximation to the entire ensemble of stereoscopically 
triangulated loops in four active regions and obtained the following results and conclusions: 

1. A linc-of-sight magnetogram that measures the longitudinal magnetic field component B z (x, y) can be 
decomposed into N m w 100 unipolar magnetic charges, from which maps of all three magnetic field 
components B x (x, y), B y (x, y), and B z (x, y) can be reconstructed at the solar surface with an accuracy 
of w 1% — 2%. This boundary condition allows us to compute a 3D potential field model B p (x) in a 
3D cube encompassing an active region, by superimposing the potential fields of each buried magnetic 
charge. Our algorithm takes the curvature of the solar surface into account and is accurate up to about 
a half solar radius away from disk center. 

2. Forward-fitting of the twisted flux tube model to 70-200 loops per active region improves the median 3D 
misalignment angle between the theoretical field lines and the observed stereoscopically triangulated 
loops from fj, = 19° — 46° for a potential field model to /i = 14° — 19° for the non-potential field 
model, which corresponds to a reduction of A/i = /i p — fi NP = 3°, 26°. The residual misalignment 
is commensurable with the estimated stereoscopic measurement error of [ise ~ 8° — 12°. 
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3. The application of our stereoscopy-constrained model allows us to obtain maps of the non-potential 
magnetic field components, B x (x, y), B y (x, y), B z (x, y), the force- free a-parameter a(x, y), and the cur- 
rent density j z (x, y) at the photospheric level or in an arbitrary 3D computation box. The divergence- 
freeness and force-freeness, numerically evaluated over a 3D computation box, was found to be reason- 
able well fulfilled. 

4. The statistics of parameters obtained from forward-fitting of our force-free model yields the following 
values: field line lengths L s=s 50 — 300 Mm (median 163 Mm), number of twist turns N tW i S t IS 0.25 
(median N tW ist = 0.06), nonlinear force-free a-parameter a < 15 x 10~ n cm -1 (median a«4x 10~ n 
cm -1 ), and current density \j z \ w 10~ 2 - 10 4 Mx cm' 2 s" 1 (median \j z \ w 1500 Mx cm -2 s _1 ). All 
twisted loops are found to be far below the critical value for kink instability (N tW i S t ~ 1-25 turns) and 
Joule dissipation of their currents (with a median Poynting flux of Fjj < 20 erg cm -2 s _1 ) is found be 
be far below the coronal heating requirement (Fjj ~ 10 5 — 10 erg cm -2 s _1 ). 

Where do we go from here? The two algorithms developed here provide an efficient tool to quickly 
compute a potential field and a quasi-force-free solution of an active region, yielding also accurate mea- 
surements of geometric (helically twisted) loop parameters (full length of field line, number of twist turns), 
based on the constraints of stereoscopically triangulated loops. Our approach of forward-fitting a magnetic 
field model to coronal structures has also the crucial advantage to bypass the non-force-free zones in the 
lower chromosphere, which plague standard NLFFF extrapolation algorithms. A next desirable project 
would be to generalize the non-potential field forward-fitting algorithm to 2D projections [x(s),y(s)] of 3D 
loop coordinates [x(s),y(s),z(s)] (as they are obtained from stereoscopic triangulation), so that the model 
works for a single line-of-sight magnetogram combined with a suitable single-spacecraft EUV image, without 
requiring stereoscopy. The stereoscopic 3D loop measurements (used here) as well as 3D vector magneto- 
graph data (once available from HMI/SDO), however, represent important test data to validate any of these 
non-potential field models. 



APPENDIX A: Deconvolution of Magnetic Charges 

A 3D parameterization of a line-of-sight magnetogram B z {x,y) can be obtained by a superposition of 
buried magnetic point charges, which produce a surface magnetic field B = (B x , B y , B z ) that is constrained 
by the observed magnetogram, as defined in Section 2.1. Here we describe the geometric inversion of the 3D 
coordinates (x m ,y m , z m ) and surface field strength B m of a single magnetic point charge M that produces 
a local peak B z with width w at an observed position (x p ,y p ) (Fig. 14, left). The geometric relationships 
can be derived most simply in a plane that intersects the point souce P and the line-of-sight axis. In the 
plane-of-sky we define a coordinate axis p that is orthogonal to the line-of-sight axis z and is rotated by an 
angle 7 with respect to the x-axis, so that we have the transformation, 

Pp = \jx 2 P +y 2 P 

x p = p p cos( 7 ) 
y p = pp sin(7 ) 

z p = \j l -Pl 

Thus we have the four obscrvables (B z ,p p , z p , w) and want to derive the model parameters (B m ,x m , y m , z m ). 
In Fig. 14 (right hand side) we show the geometric definitions of the depth d m of the magnetic charge (the 
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radial distance between M and the surface) , the distance d between the magnetic charge M and the surface 
field at an observed position P, which is inclined by an angle j3 to the vertical direction above the magnetic 
charge M, so we have the relation, 

cos/3=^. (A2) 

The line-of-sight component B z at point P has an angle (a — /?) to the radial direction B r with field strength 
B r = B m (d m /d) 2 (Eq. 1), and thus obeys the following dependence on the aspect angle a and inclination 
angle f3 (with A2), 

B z = B m (^j cos {a-P) = B m cos 2 (/3) cos (a - 0) . (A3) 
The radial coordinate p p is related to the radial coordinate p m of the magnetic charge by, 

, . , n . , sin (a — 0) , . A . 

P P = p m + dsm(a- P) = p m + d m — . (AA) 

cos p 

Further we have the geometric relationships for the aspect angle a, the 3D distance r m from Sun center, and 
the depth d m , 

a= ar ctan (p m / z m ) 

r m = VP™ + z li ■ (-45) 

dm (1 Tm) 

The observed line-of-sight component B z has a dependence on the inclination angle (3 (Eq. A3), and thus 
we need to compute the optimum angle f3 p where the component B z has a maximum, because we can 
only measure the locations of local peaks in magnctograms B z (x,y). We obtain this optimum angle [3 p by 
calculating the derivative dB z /d(3 from Eq. (A3) and setting the derivative to zero at the local maximum, 
i.e., dB z /d(3 = at /3 = (3 p , which yields a quadratic equation for tan ((3 p ) that has the analytical solution, 

tan(/3 p) = ^±^g^^. (A6) 
yHp ' 4tana 3 v ; 

For the special case of a source at disk center (a = 0) this optimum angle is (3 p = 0° , but increases monotoni- 
cally with the radial distance from disk center and reaches a maximum value of (3 P = arctan (l/\/2) w 35.26° 
at the limb (a = n/2). 

Further we need to quantify the half width w of the radial magnetic field profile B z (p) across a local 
peak, which is also one of the observables. A simple way is to approximate the magnetic field profile B z (p) 
with a Gaussian function, which drops to the half value B z /2 at an angle /3 = fc, 

B z {$ = /3 2 ) = \b x {0 = p p ) . (A7) 

Inserting the function B z ((3) (Eq. A3) yields then the following relationship for the half width w = (p2 — p p ), 

sin (a — P2) sin (a — j3 p ) 



w = (p 2 - p P ) = d. n 



cos P2 cos j3 p 



{AS) 



For the special case at disk center a = 0, the solution is w — w(ao) — d m tan (/3 2 ), and /3 2 is related to /3 P 

by 

cos3 (/32 ).£^M. (A9) 
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A general analytical solution is not possible, but numerical inversions of Eq. (A3) give the following very 
close approximation, 

w w d m tan/3 2 cos a (1 — 0.1a) . (^410) 

We have now all geometric relationships to invert the theoretical parameters {B m ,x mi y m , z m ) from 
the observables (B z , pp, zp,w). An explicit derivation of the theoretical parameters is not feasible, but an 
efficient way is to start with an approximate value for the aspect angle, since tan(a) = (p m /z m ) w (p p /z p ) 7 
followed by a few iterations to obtain the accurate value. The inversion can be done in the following order 
(using Eqs. A1-A10), 

« arctan(/9p/z p ) 



a 

P P 
B m 

$2 
dm 

Pm 

Zm 
%m 



= arctan 



(\/9 + 8tan 2 a - 3^ /4 tan a 
B z /[cos 2 P P cos(a-/3 p )] 
arccos ((cos (3 p ) 3 /2) 1/3 
wj [tan/?2 cos a (1 — 0.1a)] 
(1 - d m ) 

p p - d m s'm (a - f3 p )/ cos f3 p 

, / r 2 _ „2 
V m rm 

p m COS 7 

p m sin 7 



(All) 



In our algorithm we iteratively determine the local peaks of the line-of-sight components B z (x, y) and their 
widths w, which are found to be close to Gaussian 2D distribution functions, and invert the model parameters 
(B m , x m , y m , z m ) with the inversion procedure given in Eq. (All). 
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Table 1. Data selection of four Active Regions observed with STEREO/EUVI and SOHO/MDI. 



# 


Active 


Observing 


Observing 


Spacecraft 


Number 


Magnetic 


Magnetic 




Region 


date 


times 


separation 


of EUVI 


field strength 


flux 




Active 




(UT) 


angle (deg) 


loops 


B(G) 


(10 22 Mx) 


A 


10953 (S05E20) 


2007-Apr-30 


23:00-23:20 


6.0° 


200 


[-3134,+1425] 


8.7 


B 


10955 (S09E24) 


2007-May-9 


20:30-20:50 


7.1° 


70 


[-2396,+1926] 


1.6 


C 


10953 (N03W03) 


2007-May-19 


12:40-13:00 


8.6° 


100 


[-2056,+2307] 


4.0 


D 


10978 (S09E06) 


2007-Dcc-ll 


16:30-16:50 


42.7° 


87 


[-2270,+2037] 


4.8 



Table 2. Misalignment statistics of four analyzed active regions. 



Parameter 


2007-Apr-30 


2007-May-9 


2007-May-19 


2007-Dec-ll 


Misalignment NLFFF 1 


24° - 44° 








Misalignment PFSS 2 


25° ± 8° 


19° ± 6° 


36° ± 13° 


32° ± 10° 


Potential Field 3 


29° ±22° 


19° ± 7° 


46° ± 23° 


21° ±9° 


Force-Free Field 4 


17° ± 7° 


16° ± 6° 


19° ± 7° 


14° ± 6° 


Stereoscopy error SE 5 


9° 


8° 


12° 


9° 



1 ) NLFFF = nonlinear force- free field code (DeRosa et al. 2009), 
2 ) PFSS = Potential field source surface code (Sandman et al. 2009), 

3 ) Potential field model (this work), 

4 ) Force- free field model (this work). 

5 ) SE: Measured from inconsistency between adjacent loops. 
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Fig. 1. — Flow chart of magnetic field modeling algorithm: Starting from a line-of-sight magnetogram (top 
left), the 3D magnetic field components and field lines of a potential field are calculated. On the other side, 
3D loop coordinates are calculated from a pair of STEREO A and B images using stereoscopic triangulation 
(top right). A quasi- force- free magnetic field model with variable parameters that represent twisted loops is 
then forward-fitted to the loop coordinates by minimizing the misalignment angles between the theoretical 
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Fig. 2. — The field line geometry is shown for an untwisted cylindrical flux tube (left), a twisted cylindrical 
flux tube (middle), and for a twisted radial field (right), from the side view in the xz-plane (top) and from 
the top view in the a;y-plane (bottom). The top panels show the longitudinal magnetic field component 
B a (p) and the bottom panels show the azimuthal magnetic field component B v (p, ip) . 



-0.2 -0.1 0.0 0.1 0.2-0.2 -0.1 0.0 0.1 0.2 
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-0.2 -0.1 0.0 0.1 0.2-0.2 -0.1 0.0 0.1 0.2 



Fig. 3. — Simulations of potential field(aj = 0; left side) and twisted non-potential field cases (atj ^ 0; right 
side), represented by line-of-sight magnetograms B z (x, y) (grey scale) and magnetic field lines B(s) projected 
into the x — y plane. Each set includes a unipolar magnetic charge (first row), a dipolar (second row), and 
a quadrupolar configuration (third row). 
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-0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 



MDI20070519J240 Model map (N m =200) Difference 
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MDI20071211J630 Model map (N m =200) Difference 
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Fig. 4. — Left column: photospheric magnetograms observed with SOHO/MDl for four observations. Middle 
column: the magnetic held is decomposed into N m = 200 unipolar magnetic charges and the model map 
displays the line-of-sight component of the magnetic held, B z . Right column: difference between the observed 
MDI magnetogram and the model. 
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peak(P) =29.1+21.6 
median (M) = 29.7 
average (A) = 39.9 



20 40 60 80 

Misalignment angle a mis (deg) 



Fig. 5. — A number of 200 stereoscopically triangulated loops observed with STEREO (blue) are compared 
with the theoretical field lines of a potential field (red curves), defined by the potential fields of N m = 100 
buried magnetic charges, decomposed from the MDI line-of-sight magnetogram. The 3D views are shown in 
a top-down view in the x — y plane (middle), and in sideviews in z — y plane (right) and x — z plane (top). 
A histogram of the 3D misalignment angles is shown in the bottom panel, with the average (A), median 
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MDI20070430 2300Aa1 00 : div-free=0.000050, force-free=0.000536 
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ngrid = 1 
thresh = 0.0 G 
imag/nmag=100/100 
nsm = 3, iopt = 
qloop = 0.50 
nseg = 10 
acc = 0.0010 
amax = 100.0 
iter/niter= 1 5/ 20 
qzone = 0.500 
qv = 0.0 
cpu = 1257.4 s 




0.8 0.9 1.0 1.1 



peak (P) =16.7+ 6.7 
median (M) = 19.7 
average (A) = 29.7 



20 40 60 

Misalignment angle a mis (deg) 



Fig. 6. — Similar representation as in Fig. 5 for active region A (2007 Apr 30), but modeled by forward- 
fitting our analytical approximation of a force-free field model. Note that the 3D misalignment angles for 
this quasi-force-free field are significantly smaller than for the potential-field model shown in Fig. 5. See also 
movie A on forward-fitting of active region A. 
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Fig. 7. — Similar representation as in Fig. 6, for active region NOAA 10955 observed on 2007 May 9, 20:30 
UT. See also movie B on forward-fitting of active region B. 
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MDI20070519_1240Aa100 : div-free=0. 143793, force-free=0. 187655 
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40 60 
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Fig. 8. — Similar representation as in Fig. 6-7 for active region NOAA 10953 observed on 2007 May 19, 12:40 
UT. See also movie C on forward-fitting of active region C. 



-32 - 



MDI20071211 1630Aa100 : div-free=0.000095, force-free=0.000134 
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40 60 
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Fig. 9. — Similar representation as in Fig. 6-7 for active region NOAA 10978 observed on 2007 Dec 11, 16:30 
UT. See also movie D on forward-fitting of active region D. 
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Fig. 10. — The distribution of best-fit misalignment angles a mis (between the observed loops and the best 
fit of the theoretical force-free field model) is shown (left panels). The distributions are characterized by the 
value of the peak (P) of a Gaussian fit with Gaussian width, by the median (M), and the average (A). The 
dependence of the average misalignment angle (A) as a function of the number of magnetic components N m 
is shown in the right panels, for the potential field model (P) and the non-potential (force-free field) model 
(NP). 




Fig. 11. — The magnetic field component maps B p (x, y), B p (x, y), and B P (x, y) are shown for the potential 
field model (left column), difference maps of the non-potential field model to the potential field model, 
B^ p — B p , By P — B p , and B^ p — B P (middle column), the observed line-of-sight component B° bs (x,y) 
(bottom right panel), and the derived nonlinear-force free a-parameter map a(x,y) (top right panel), and 
electric current j z (x, y) map (middle right panel), for active region A (2007 Apr 30). The grey scale is identical 
for the potential field and difference maps. The number of magnetic charge components is N m — 100. 




Fig. 12. — Observed SOHO/MDI magnetogram B z (x, y) (top row) and electric current maps j(x, y) obtained 
from the force-free field forward-fitting, for active regions B (2007 May 9), C (2007 May 19), and D (2007 
Dec 11). 
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Fig. 13. — Statistical distributions of loop-associated field line lengths L (top row), the number of twist 
turns per loop, N tw i s t (second row), the force- free a-parametcr (third row), and the current density \j z \ 
(bottom row), in lin-lin (left column) and log- log histograms (right column), with the median indicated 
(dashed vertical lines) and powerlaw fits of the distributions. 
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Fig. 14. — 3D geometry of a point source P = (x p ,y p ,z p ) in a cartesian coordinate system is shown (left), 
with the z-axis aligned to the line-of-sight from Earth to Sun center. The plane through the line-of-sight 
axis and the point source P has a position angle 7 in the plane-of-sky with respect to the x-axis and defines 
the direction of the axis p. The geometry of a line-of-sight magnetic field component B z is shown in the 
(z, p)-plane on the right hand side. A magnetic point charge M is buried at position (z m ,p m ) and has an 
aspect angle a to the line-of-sight. The radial component B r is observed on the solar surface at location P 
and has an inclination angle of f3 to the local vertical above the magnetic point charge M. The line-of-sight 
component B z of the magnetic field has an angle (a — (5) to the radial magnetic field component B r . See 
details in Appendix A. 



